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Abstract 


The  paper  is  an  invited  presentation  at  the  Sixth  IMAC  International 
Symposium  on  Computer  Methods  for  Partial  Differential  Equations,  June  23-26, 
1987,  Lehigh  University.  It  discusses  computational  complexity  of  various 
versions  of  the  finite  element  method  in  relation  to  the  achieved  accuracy  of 
the  finite  element  solution. 
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1.  Introduction 

There  are  many  veraiona  and  forma  of  finite  element 
methoda  implemented  in  the  finite  element  codea.  The 
baaic  architecture  of  the  etandard  finite  element  pro- 
grama  are  eaaentially  very  aimilar.  Once  the  claaa  of 
problema  that  can  be  aolved  haa  been  aelected  (elaatic- 
•*y.  plates  and  shell  problems,  heat  transfer  etc.),  one 
of  the  main  featurea  ia  the  variational  form  given  to  the 
problem  (primary,  mixed,  hybrid  etc.)  together  with  the 
type  of  finite  element  used  (conforming,  non-conforming 
etc.). 

Nowadays  several  hundred  of  finite  element  program 
ay  sterna  are  available,  covering  a  wide  range  of  claaa 
of  problems.  Despite  the  great  number  of  codes  it  ia 
very  difficult  to  find  comparison  results  between  differ¬ 
ent  programa.  The  main  reason  is  due  to  the  fact  that 
the  comparison  is  a  very  complex  task  involving  a  large 
number  of  different  factors.  The  usual  approach  ia  to 
test  different  methoda  and  codes  on  benchmark  prob¬ 
lems  (see,  e.g.,  MacNtal  and  Harder  1 1 4 j ,  Robinson 
and  Blackham  (20,2 1 1)  -  Although  this  approach  has  ob¬ 
viously  various  drawbacks,  it  gives  good  and  concrete 
data  for  an  effective  comparison.  Of  course  additional 
aspects  like  convergence,  rate  of  convergence,  effect  of 
the  computer  technology  etc.,  are  essential  as  well.  Dif¬ 
ferent  simple  model  problems  together  with  a  theoretical 
modeling  are  also  very  worthwhile. 

The  goal  of  this  note  is  to  address  some  results  al¬ 
lowing  basic  essential  comparison  between  finite  element 
methods.  There  are  three  major  version  of  the  finite  ele¬ 
ment  method:  the  standard  h  -  version,  the  p  -  vtreton 
and  the  h-p  vtreton.  The  p-  and  h-p  versions  are 
a  recent  development.  In  Ssabo'  [2S|,  Guo  and  Babuika 
|10|,  Seri  [22)  recent  advances  in  the  p-  and  h-p  ver¬ 
sions  are  given.  For  a  survey  of  the  state  of  the  art 
see  also  Babulka  (3 [ .  We  will  try  to  make  some  compar¬ 
isons  between  these  three  versions,  focusing  especially  on 
their  computational  aspects.  We  will  particularly  ana¬ 
lyse  the  relations  between  the  accuracy  of  the  error  and 
the  computational  work  Of  course  any  comparison  is 
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very  complex,  hence  we  will  consider  only  some  aspects 
of  finite  element  method  for  linear  problema.  Finally  we 
will  consider  a  benchmark  problem  and  we  will  see  that 
the  practical  results  are  in  good  agreement  with  the  ones 
based  on  the  computational  theoretical  models. 

X.  The  scheme  of  the  finite  element  method 

The  finite  element  method  consists  of  few  bask 

phases: 

(a)  Topol  op  y 

By  this  we  mean  the  mesh  generation.  The  mesh 
has  to  respect  the  geometry  of  the  domain  under  con¬ 
sideration  and  the  requirement  that  effective  and  accu¬ 
rate  solution  will  be  obtained.  The  mesh  generation  is 
a  very  laborious  part  of  the  method  also  when  sophis¬ 
ticated  mash  generators  are  available.  Very  often,  es¬ 
pecially  in  three  dimensions,  the  mesh  reflects  only  the 
geometry. 

* 

(b)  Local  etif  fntee  matrtcee 

The  complexity  of  the  computation  of  local  stiffness 
matrices  depends  on  the  degree  of  the  elements  and  the 
type  of  hardware  used  (sequential  or  parallel). 

Let  us  consider  here  computations  of  a  C“  quadrilat¬ 
eral  element  of  degree  p  for  a  two  fields  (elasticity  in  two 
dimensions)  problem.  Using  hierarchical  elements  (see. 
e.g.,  Noor  and  0ak«Jfcs|4|,  Ssabo'  |23|),  an  element  of 
degree  p  has  Q  shape  functions,  where 

(1)  <3  =  8  .  P  =  I 

Q  =  2  I  4p  +  max(0,  j(p  2)(p  -  3))  .  p  >  I 

The  second  part  in  the  formula  for  Q  gives  the  number 
of  interns/  shape  functions  which  are  treated  later  by 
condensation  procedure. 

Proper  programmming  on  a  sequential  machine  can 
be  made  so  that  0(p*)  operations  are  needed  for  com¬ 
putation.  An  experimenetal  program  in  the  (normalised 
VAX)  time  units  leads  to  the  computational  work  for  the 
local  stiffness  matrices  WLM  ,  where 


(3) 


WL 


—  (2.5  +  0.032  p4). 
25 


We  recall  that  in  th«  case  of  uniform  square  mesh 
the  number  N  of  degrees  of  freedom  is  the  following: 


Figure  1  shows,  for  different  values  of  p,  the  time  units 
needed  for  the  computation  of  25  stiffness  matrices. 
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(4)  N  =  8(m  +  l)a  ,  p=l. 

N  =2((m  +  l)*  +  2m(m  +  l)(p  -  1) 

+  m*  max(0,  i(p  -  2)(p  -  3))|,p  >  1. 

(d)  Postprocessing 

By  this  ere  mean  computations  of  desired  values, 
srror  control,  graphics  etc.  This  phase  is  often  computa¬ 
tionally  very  expensive.  Because  the  coot  of  this  phase  is 
nearly  the  same  for  all  the  versions  of  the  finite  element 
method,  we  will  not  consider  this  phase. 

We  consider  as  total  work  of  computation  W  only 
the  work  needed  for  the  phases  (b)  and  (c).  We  have 

(5)  W  m  m’  WLM  +  W ,  . 


(c)  AssemMg  end  elimination 

The  work  needed  for  this  phase  depends  on  the 
topology  and  the  degree  of  the  elements.  The  practi¬ 
cally  moat  complex  topology  in  two  dimensions  is  the 
case  of  the  uniform  square  mesh.  Assuming  a  decompo¬ 
sition  with  m*  elements,  the  experimental  program  in  the 
(normalised  VAX)  time  units  of  the  elimination  methods 
leads  to  the  following  expression  for  the  computational 
work  Wm  (in  the  range  1  <  p  <  8,  1  <  m  <  5): 

(3)  Wm  =  2  +  (0.013«  +  0.004  m’)  p’  . 

Figure  2  shows  the  experimental  time  (in  VAX  units)  as 
a  function  of  p  and  m. 


Sometimes  different  methods  are  compared  with  re¬ 
spect  to  the  number  N  of  degrees  of  freedom  instead  of 
computational  work.  Hence  we  define  a  correction  coef¬ 
ficient 
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3.  The  goal  of  the  computations  and  the  princi¬ 
ples  of  comparisons  of  the  method 


As  we  said  before,  comparison  of  methods  and  codes 
is  a  very  complex  problem.  It  depends  on  many  factors, 
one  important  being  the  man  power  cost  needed  to  op¬ 
erate  the  code  in  practical  environment  The  goal  of  the 
computation  is  to  achieve  the  desired  data  in  the  range 
of  a  given  accuracy.  We  will  try  to  assess  the  various 
versions  of  the  finite  element  method  and  we  will  con¬ 
centrate  only  on  the  energy  norm  error  c  ,  To  get 
reasonably  valid  conclusions  we  will  use 


Fig  2 


a-  theoretical  model  analysis; 
b-  benchmark  example  analysis. 

The  main  idea  of  the  theoretical  model  analysis  is 
to  assume  that  the  energy  norm  error  can  be  given  as  an 
explicit  function  depending  on  both  number  and  degree 
of  elements.  We  suppose  that 

(•)  lUlls  *  *(m,p)  , 

where  *(m, p)  is  a  known  function.  We  will  derive  suit¬ 
able  express  ions  for  the  function  ¥(m,p)  so  that  it  has 
the  forms  close  as  much  as  possible  to  the  theoretical  es¬ 
timates  for  various  versions  of  the  method.  Then  we  will 
study  the  dependence  between  the  erTor  ||e||*  and  the 
computational  work  W .  Although  this  approach  has,  of 
course,  various  shortcomings,  nevertheless  it  can  lead  to 
valuable  conclusions. 

The  benchmark  example  will  be  the  'classical*  en¬ 
gineering  problem  of  a  simply  supported  rhombic  plate. 
We  selected  to  present  only  this  benchmark  problem  be¬ 
cause  it  has  different  character  (C‘  elements  are  used) 
and  our  conclusions  based  on  plane  elasticity  model  com¬ 
putations  (C°  elements)  are  still  valid.  Of  course  ad¬ 
ditional  informations  and  estimates  could  be  here  used 
too.  Nevertheless  we  will  see  that  the  results  are  in  good 
agreement  with  the  conclusions  based  upon  the  theoret¬ 
ical  model  for  the  plane  elasticity. 

4.  Complexity  of  the  finite  element  computation 

The  effectiveness  of  the  method  depends  on  the  claes 
of  problem  characterising  the  function  ♦  (*»,*)  in  (8), 
which  is  a  form  of  an  estimate. 

a)  Tk*  |««n  uniform  mttk 

Let  u0  be  the  exact  solution  of  the  given  problem 
and  e  be  the  error  in  the  h  -  p  version.  Then  in  BabuMka 
and  5uri  |6|  the  following  estimate  has  been  proved: 

(7)  11*11#  —  pj,  _  ,  ^(*)  1 1 **o  1 1  w  *  i  o  i  ■ 

where  h  is  the  mesh  sise,  m  =  min(p,k  -  1),  fl  is  the 
given  domain,  ff*(0)  the  usual  Sobolev  space  and  C(k) 
a  constant  independent  of  h  and  p  but  depending  on  k 
and  fl.  In  our  case  we  assume  that  h  =  m'  1 

Usually  the  function  w,  can  be  imbedded  in  various 
Sobolev  spaces.  Hence  we  will  assume  that 

(•)  C(k)  ||«o  l[*  ,  —  ♦  (*)  , 

where  4(fc)  is  an  a-priori  known  function  We  can  as¬ 
sume  that 


This  type  of  function  *(fc)  expresses  the  case  when  the 
function  s  is  uniformly  untmootk  over  the  entire  do¬ 
main  or  has  singular  behaviour  in  an  a-priori  unknown 
location. 

Another  form  we  will  consider  is  the  class  of  analytic 
solutions.  In  this  case  we  will  assume 

(10)  *(*)  =  k!  f  , 

where  d  >  1  (the  coefficient  d  expresses  the  site  of  the 
natural  domain  of  the  analyticity). 

Assuming  the  computational  work  W  given  by  (5) 
and  the  function  #(k)  by  (0),  we  show  in  the  figures 
3,4 ,5,8  the  dependence  of  the  accuracy  of  the  error  ||e||* 
on  the  computational  work  W  (in  log  -  log  scale)  for 
different  values  of  ka  and  a  =  1.  The  values  shown  are 
the  minimal  errors  for  integers  k  -  1,2,., .,  k0  -  1  (of 
course  also  the  error  in  the  entire  range  for  k  could  be 
easily  computed,  but  it  would  not  change  the  conclu¬ 
sions).  In  all  these  figures  the  lines  corresponding  to  the 
values  p  *  1,2, ...,  8  are  reported.  On  each  line  different 
values  of  m  =  1,2, ...  are  considered. 
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Let  oa  now  consider  a  class  of  problems  where  the 
eolation  is  singular  in  a  vertex  of  the  domain.  In  this 
caae  the  solution  has  the  leading  singularity  of  the  type 

(n)  Is  =  '*  ♦(*)  . 

where  r  >  0,  a  >  0  and  +  is  an  analytic  vector  function 
in  #.  Then  in  Bmkutka  and  Sun  [6]  it  has  been  proved 
that 


("»"*  . 
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Let  us  now  consider  the  case  when  #(k}  is  given  by 
(10).  The  figures  7  and  S  show  the  behaviour  of  the  error 
against  the  computational  work  for  two  different  values 
of  d. 
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The  figures  0,10,1 1  show  the  accuracy  of  the  error  against 
computational  work  for  different  values  of  a. 
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So  far  we  have  mentioned  only  decomposition*  of 
the  given  domain  with  uniform  me*h.  If  the  (olution  i* 
of  the  type  (II)  and  an  optimal  mesh  is  used,  then  the 
error  behaviour  is  close  to  the  case  of  the  smooth  solution 
previously  described. 

b)  The  strongly  refined  mesh  for  the  h  -  p  version 

Here  we  consider  the  solution  u0  belonging  to  the 
space  BJ( fl)  (see  Guo  and  Babvika  [11,12]).  In  this 
case  the  following  estimate  holds  for  the  error: 

(13)  ||e||,  =  C  exp(—b  Af  • )  . 

The  topology  of  the  (strongly  refined)  mesh  is  equivalent 
to  a  rectangle  (not  a  square).  Nevertheless  we  will  use 
our  result  for  the  square  mesh  topology  (which  leads  to 
the  rate  exp(-4  S  •)).  In  this  case  the  degree  p  of  the 
elements  that  are  used  depends  on  m.  We  will  use  for 
the  error  the  relation  which  leads  to  the  above  estimate 
and  is  used  in  (11,12): 


where  [«m]  denotes  the  integer  part  of  sm.  The  figure  12 
shows  the  error  against  the  work  when  the  relation  (14) 
is  used  with  a  =  0.5. 
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c)  Conclusions 

A  careful  analysis  of  the  figures  previously  shown 
leads  to  the  following  conclusions: 

1-  Only  for  very  low  required  accuracy  or  very  un¬ 
smooth  solutions  ( k0 ,  a  small)  the  lower  order  elements 
are  preferable,  on  a  uniform  mesh.  However,  this  is  not 
usually  the  practical  case. 

2-  The  increase  of  the  quality  of  the  solution  from 
p  =  1  to  p  =  2  is  very  significant  and  p  =  2  gives  some¬ 
times  reasonable  errors  for  low  but  already  acceptable 
engineering  accuracy. 

3-  It  is  safer  to  use  elements  of  higher  degree  (3-5) 
then  lower  ones  for  the  effective  computation.  The  use  of 
higher  degree  elements  in  general  is  more  robust.  Let  us 
remark  that  in  the  elasticity  problem  the  locking  effect 
(for  v  —  0.5)  is  completely  eliminated  when  p  >  4  (see 
5cott  and  V ogelius  [27]).  This  is  another  reason  to  use 
p  >  3. 

4-  If  the  solution  is  smooth  or  the  solution  has  sin¬ 
gularities  but  the  mesh  is  properly  refined,  higher  order 
elements  are  preferable  except  for  very  low  desired  accu¬ 
racy. 

5-  The  most  effective  way  is  to  combine  the  degree  of 
the  elements  with  a  properly  designed  mesh.  Of  course 
the  effectiveness  depends  very  much  on  the  proper  design 
of  the  mesh.  For  questions  of  this  type  we  refer  to  Guo 
and  Babutka  [13],  Rank  and  Babvika  [18]. 

6-  A  high  practical  accuracy  can  usually  be  obtained 
with  the  elements  of  degree  p  <  8  when  reasonable  mesh 
is  used  in  all  practical  cases. 


(14)  11*11.  =  exp(-ap)  ,  p  =  [«m|  , 


S.  A  benchmark  problem:  the  rhombic  plate 

As  we  said  in  the  introduction,  the  usual  way  to 
compare  different  finite  element  methods  is  to  test  them 
on  a  benchmark  problem,  where  the  solution  is  in  some 
way  explicitly  known.  In  order  to  evaluate  the  effective¬ 
ness  of  the  theoretical  error  analysis  model  previously 
introduced  several  teats  has  been  performed.  In  particu¬ 
lar  we  refer  to  BabuSka  and  Suri  [6]  for  various  results 
related  to  the  solution  of  the  elasticity  problem  of  an 
L-shaped  domain.  However,  we  show  that  our  conclu¬ 
sions  are  useful  for  more  general  situations  beyond  the 
specific  one  from  which  it  has  been  derived  (of  course  a 
similar  model  could  be  made  based  on  available  results 
for  plate  bending  problems).  Hence  we  have  analyzed 
another  "classical”  benchmark  problem,  the  bending  of 
a  simply  supported  rhombic  plate  under  normal  loading, 
to  evaluate  the  generality  of  our  approach. 

a)  The  rhombic  (Kirehhof  f)  plate  problem 

The  most  severe  test  for  plate  elements  is  the  so- 
called  Morley’s  skew  plate  (see  Morley  [15,16]),  which 
is  a  uniformly  loaded  and  simply  supported  plate  of  the 
shape  of  an  equilateral  parallelogram  (rhombus) . 

In  Fig. 13  a  rhombic  plate  is  shown  together  with  a  2  x  2 
mesh  for  triangular  elements. 


Fig. 13 


Because  of  the  stress  singularity  at  the  obtuse  cor¬ 
ners  of  tha  plate  the  convergence  of  the  finite  element 
solutions  has  been  found  to  be  extremely  slow.  The  sin¬ 
gularity  depends  on  the  angle  6  and  as  6  decreases  the 
singularity  becomes  stronger.  Robinson  [19]  has  com¬ 
piled  the  results  of  many  finite  element  codes  for  this 
problem.  The  results  of  most  of  the  elemnts  show  a  very 
large  error  for  the  center  deflection  of  the  plate  even  when 
very  fine  mesh  is  used.  However,  no  indications  are  given 
about  the  real  cost  of  using  each  element  to  get  a  desired 
accuracy.  Our  main  goal  is  to  test  some  elements  on  the 
rhombic  plate  taking  into  account  accuracy  and  compu¬ 
tational  work. 


(b)  The  finite  elements 

Several  finite  elements  for  plate  are  known.  We  have 
restricted  our  attention  to  conforming  elements,  that  is 
elements  satysfying  the  Cl  -condition,  allowing  the  dis¬ 
crete  approximation  space  to  be  included  in  the  contin¬ 
uous  one.  The  main  reason  is  that  conforming  elements 
have  some  good  properties  of  monotonicity  (i.e.  for  the 
energy)  enabling  an  effective  control  of  the  error.However 
in  many  codes  non-conforming  elements  are  used. 

The  elements  we  have  tested  are  the  following: 

1-  reduced  Hsieh-Clough-Tocher  triangle  (HCTR), 

2-  assumed  stress  hybrid  triangle  (HYBR), 

3-  Argyris  triangle  (ARGY). 

The  Fig. 14  shows,  with  the  usual  convention  for  no¬ 
tation,  the  type  of  degrees  of  freedom  associated  with 
each  type  of  element.  For  details  we  refer  to  [8]  for 
HCTR,  [7,17]  for  HYBR,  [1]  for  ARGY. 
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Fig. 14 


We  only  recall  that  the  following  abstract  error  es¬ 
timates  hold  provided  a  sufficiently  smoothness  of  the 
solution: 

ll«!ls  <  C  h  ||u||s 
for  both  HCTR  and  HYBR  elements, 

ll«ll«  <  C  h*  ||u||. 

for  ARGY  element. 

(c)  The  numerical  results 

Numerous  tests  have  been  performed  for  different 
shapes  of  the  plate  and  for  different  boundary  conditions 
Here  we  focus  our  attention  on  the  energy  error  obtained 
for  various  values  of  the  angle  6.  We  have  considered  a 
plate  with  the  following  data: 


•id*  length  a  =  1.0  in, 
thickness  t  =  0.01  in. 

Young’*  modulus  E  =  30  x  10*  16/m , 

Poisson’s  ratio  v  —  0.3 
normsl  load  p  =  1.0  lb/ in'1 . 

For  the  piste  problem  we  have  computed  the  (rela¬ 
tive)  energy  error,  defined  in  the  following  way: 

ll'ila  -  ( *•-='  ~)  f  >■  >M  ■ 

where  EM,  Ec  represent  the  exact  and  the  computed 
energy,  respectively.  The  computation  has  been  carried 
out  using  routines  of  the  code  MODULGF  (see  [9])  and 
the  final  linear  system  has  been  solved  by  Cholesky  pro¬ 
cedure.  The  computation  was  performed  on  an  Apollo 
DN300  computer.  We  have  used  existing  subroutines  to 
avoid  any  bias  on  the  reached  conclusions. 

The  Fig. 15, 16  show  the  behaviour  of  the  energy  er¬ 
ror  against  the  total  time  of  computation  for  different 
values  of  the  angle  S.  As  expected,  the  error  increases  as 
the  angle  6  decreases.  However  we  note  that  in  all  the 
situations  the  higher  degree  element  ARGY  gives  bet¬ 
ter  performances  over  the  lower  degree  elements  HCTR 
and  ARGY.  More  important,  if  we  compare  these  results 
with  the  Fig.9, 10,11  we  see  the  agreement  between  our 
model  and  the  concrete  numerical  results. 

To  better  understand  the  equivalence  we  show  in 
Fig. 17, 18  the  lines  corresponding  to  degrees  p=  1  (HCTR 
and  HYBR  element)  and  p=t  (ARGY  element)  for  values 
of  <*=0.8  and  a=0.2  based  on  computation  using  relation 
(12)  (in  the  same  way  as  Fig.9,10,11  were  computed). 
These  values  of  the  parameter  describing  the  singular¬ 
ity  are  theoretically  related  to  the  values  of  6=80  and 
6=30,  respectively.  The  direct  comparison  (Fig. 15  and 
17;  Fig. 16  and  18)  shows  a  good  agreement  between  our 
model  and  practical  results. 
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0.  Conclusion! 

1-  Models  of  computational  complexity  based  on 
theoretical  error  estimates  for  typical  classes  of  solutions 
and  measured  computer  timing  of  benchmark  problems 
leads  to  reliable  assesment  and  conclusions  for  computa¬ 
tional  comparisons  of  the  finite  element  methods. 

2-  Higher  order  elements  are  preferable  because  they 
usually  need  smaller  computational  effort  for  the  same 
accuracy  of  error.  Moreover,  high  degree  elements  are 
much  more  robust  (e.g.  also  with  respect  to  the  locking 
problem).  The  less  favourable  effectiveness  in  compari¬ 
son  with  lower  order  elements  in  extreme  situations  is 
significantly  outweighed  by  very  good  performance  in 
normal  situations,  especially  using  a  "reasonable”  mesh. 

3-  The  concrete  conclusions  previously  listed  down 
in  section  4c)  are  valid  for  large  classes  of  problems. 

4-  Our  comparisons  were  based  on  the  performance 
with  respect  to  the  energy  norm.  We  expect  the  same 
conclusions  for  other  norms  and  performance  measures 
although  a  larger  emphasis  could  be  necessary  on  the 
used  mesh  (e.g.  dealing  with  the  pollution  problem;  see 
BabuSka  and  Oh  [5]). 

5-  We  did  not  address  the  question  of  man  power 
cost.  For  the  comparison  of  this  type  in  industrial  envi¬ 
ronment  we  refer  to  Barnhart  and  H ciscmann  [2]  where 
higher  order  elements  were  shown  highly  preferable. 

6-  We  did  not  address  here  the  question  of  the  rela¬ 
tion  of  the  quality  results  assessment  (see,  e.g.,  Noor  and 
BabuSka  |4|)  to  the  degree  of  the  elements.  This  field  is 
quite  new  for  all  methods.  Nevertheless  hierarchical  in¬ 
crease  of  p  leads  to  a  simple  and  effective  error  control 
for  all  data  of  engineering  interest  baaed  on  various  types 
of  extrapolation  (see,  e.g.,  (23,24,28]). 

7-  We  did  not  address  the  effectiveness  of  the  finite 
element  treatment  of  various  mathematical  formulations 
of  the  same  problem  (for  example  Kirchhoff  and  Mindlin- 
Reissner  models  for  plate  etc.).  This  comparisons  will  be 
made  in  a  forthcoming  paper. 
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